Noise-induced metastability in biochemical networks 
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Intra-cellular biochemical reactions exhibit a rich dynamical phenomenology which cannot be 
explained within the framework of mean-field rate equations and additive noise. Here, we show that 
the presence of metastable states and radically different timescales are general features of a broad 
class of autocatalytic reaction networks, and that this fact may be exploited to gain analytical 
results. The latter point is demonstrated by a treatment of the paradigmatic Togashi-Kaneko 
reaction, which has resisted theoretical analysis for the last decade. 
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With recent advances in experimental techniques, it is 
becoming increasingly clear that the dynamics of cellu- 
lar biochemical reactions are subject to a great deal of 
noise [IJ . This poses a significant challenge to our under- 
standing of such systems, as it has been known for some 
time that the effects of noise may lead to substantial dif- 
ferences in the macroscopic behavior 0, [1] • The reactions 
which take place within a cell are highly interdependent, 
together forming biochemical networks which support the 
functioning of the cell. It remains a major open problem 
to make clear the link between the structural features 
of these networks and the resulting dynamics. A full 
understanding of the effects of noise is essential to this 
effort 

Here, we report analytical progress on this problem 
made by studying a simple class of autocatalytic reaction 
networks whose dynamical behavior is radically affected 
by intrinsic stochasticity in finite volume cells. In partic- 
ular, we show how networks of this type give rise to a sep- 
aration of timescales between fast almost-deterministic 
oscillations and slow stochastic metastability. Our class 
includes the influential Togashi-Kaneko (TK) reaction 
scheme, numerical simulations of which have been found 
to undergo a noise-induced dynamical transition [el. De- 
spite the importance of their work, a satisfactory analytic 
treatment of this effect has not been achieved in over a 
decade. Here we provide such a treatment as an applica- 
tion of our theory. 

The general model we work with is composed of n 
chemical species, denoted by Xi with i = re- 
siding in a cell of (non-dimensional) volume V . The 
molecules undergo autocatalytic reactions of the form 
Xi + Xj — )■ 2Xj , with rate coefficients . We put = 
if that particular reaction is not possible. We also stip- 
ulate that the total rates of creation and destruction of 



FIG. 1. (Color ordine) Sample stochastic time series of the 
simple three-species reaction network described in the text, 
with volume V — and diffusion coefficient D = 10"'*. 
The thick (blue), thin (red) and dashed (purple) lines show 
the concentrations of chemicals Xi, X2 and X3, respectively. 
The smaller figures show detail of rapid oscillations (left) and 
metastability (right), taken from the main plot. All simula- 
tions were performed using Gillespie's algorithm [^. 



components Xi — Xi/V . 

The dynamics of the system defined by the above re- 
actions are specified once the transition rates, T{x\x'), 
indicating the probability per unit of time that the sys- 
tem goes from state x' to state x, are given. They are 
found by invoking mass action: 



/ 1 1 



{ 



DVx. 



Tlx, 



1 
V 



(1) 



DV . 



The probability of finding the system in the state x at 
time t, P{x,t), then satisfies the master equation 



each reactant i are in balance, that is, r^- — ^ 
Two additional reactions, ^ Xi and Xi — ^ 0, repre- 
sent diffusion into and out of the cell, respectively. The 
rate of diffusion is slow compared to the internal reac- 
tions, having coefficient D <C 1. We will also use the 
symbol Xi to denote the number of molecules of that 
type, and x to indicate the concentration vector with 
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with the transition rates given above [7|. 

Stochastic simulations of reaction networks of the class 
described above display a rich phenomenology includ- 
ing rapid oscillations and random switching between 
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metastable states. For example, the time series dis- 
played in Fig. [1] were obtained from simulations of a 
three-species reaction with (arbitrarily chosen) non-zero 
reaction rates ri 2 — 1, ?'2,3 = 4, r3_2 = 3, r^ i = 1. In 
what follows we will show how these features can be qual- 
itatively and quantitatively understood by an analysis of 
the influence of noise and the separation of timescales. 

The dynamics are drastically affected by the relation- 
ship between the cell volume and the diffusion coeffi- 
cient. To elucidate this, we introduce a rescaled volume 
A = DV, which we treat as an 0(1) control parame- 
ter. Scaling V and D simultaneously in this way, we 
can rewrite the master equation ([2]) as a power series 
in a single small parameter (we choose D, but is 
also a valid expansion parameter) , leading to a Kramers- 
Moyal expansion 1)]. Truncating it at second order, one 
obtains a Fokker-Planck equation equivalent to the fol- 
lowing stochastic differential equation (SDE), defined in 
the Ito sense 0: 



Xi — X' 



J2R^]Xi+Dil-x,) + ^rj,{t), (3) 



where i — I, . . . ,n, Rij — rji — rij and the rji are Gaussian 
noise variables with zero mean and correlator 
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Here the angle brackets signify an average over the noise, 
and Sij = rij -\- Vji. 

Several important facts about the dynamics can be as- 
certained from inspection of Eqs. ([3]) and First, we 
discuss the limit of large volume. The factor of A~^ in 
the noise correlator indicates that for finite volumes the 
system experiences internal fluctuations. These vanish 
as A — >■ 00, leaving behind a deterministic system of dif- 
ferential equations equivalent to those obtained from a 
mean-fleld analysis of the reaction network. For general 
reaction networks these equations describe simple oscil- 
latory relaxation towards the homogeneous fixed point 
Xi = 1 for all i. This prediction is quite at odds with 
the rich phenomenology which is observed in stochastic 
simulations (as seen in Fig. [U for example). A proper 
treatment of the noise is thus necessary: from now on we 
keep A fixed and finite. 

The presence of the small parameter D in Eq. ^ im- 
plies a separation of timescales. On an 0(1) timescale 
(which we refer to as fast) , diffusion is negligible and the 
system feels no noise. Setting D = in Eq. Q yields 
a deterministic dynamical system in which the homo- 
geneous state = 1 is a center; it has Jacobian ma- 
trix R, which is antisymmetric and thus has all imagi- 
nary eigenvalues. We can therefore expect rapid almost- 
deterministic oscillations as seen, for example, in the 



lower left panel of Fig. [TJ On a slow 0{1/D) timescale, 
two additional factors play a role. First, the system expe- 
riences a deterministic linear drag towards the homoge- 
neous state. Second, the effects of noise become relevant, 
leading to stochasticity in the trajectories. 

For smaller volumes, the overall noise strength is 
greater, and thus the form of the noise correlator (|4]) 
has an important role in shaping the system dynamics. 
In particular, since the strength of the noise is a function 
of the state of the system, trajectories are forced away 
from states giving rise to large values of noise, creating 
an effective attraction towards those states in which the 
noise vanishes. This effect is relatively well-known in the 
study of systems with multiplicative noise (for example, 
see [10[ and references therein), and we will illustrate it 
with an explicit calculation for the TK model. Inspection 
of the correlator ^ reveals that the states for which the 
noise vanishes are those in which no autocatalytic reac- 
tion can occur. That is, for each pair i,j one of Xi, xj 
or rij must be zero. The metastability of these states 
is further enhanced by the fact that this condition also 
causes the 0(1) term in Eq. ^ to vanish. An example 
can be seen in the lower right panel of Fig. [l] where the 
state Xi — 3, X2 = , X3 = is metastable. 

As well as providing a qualitative picture of dynamics 
observed in this class of biochemical reaction networks, 
the mathematics we describe may also be employed to 



(4) obtain precise analytical results We now illustrate 



these methods in the paradigmatic case of the TK reac- 
tion 0. The model is composed of four chemical species 
whose reactions form a closed cycle, so that the non-zero 
rates are ri 2 = ^"2,3 = r^^^ = r^ i = 1. In stochastic sim- 
ulations of the model, different dynamics are observed 
depending on the volume of the cell. For very large vol- 
umes, one finds an approximately homogeneous distribu- 
tion of chemical species, however, at lower volumes the 
system is typically dominated by a pair of species (either 
Xi and X3, or X2 and X^), with the other pair absent: 
these are the metastable states predicted in the earlier 
discussion. 

To visualize this dynamical transition, TK ^ intro- 
duced the quantity z = {xi + x^) — {x2 + X4). The 
pair-dominated state corresponds to \z\ « 4. By measur- 
ing the stationary distribution P{z) from long simulation 
runs, one observes a transition induced by cell volume - 
see Fig. [21 There is a critical volume K ~ at which 
P{z) is flat; above Vc the distribution has a single peak 
at z = 0; below Vc it is bimodal with peaks at z w ±4, 
indicative of the pair-dominated regime. 

In large volumes the model also exhibits quasi-cycles, a 
second (weaker) stochastic effect whereby damped oscil- 
lations present in the deterministic dynamics are excited 
by the noise. Quasi-cycles are amenable to analysis using 
a linear noise approximation [12|, however, it is clear that 
the dynamical transition is related to the noise-induced 
metastability discussed above and will require more pow- 
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FIG. 2. (Color online) Stationary probability distribution for 
z = {xi + X3) — {x2 + Xi) in the TK reaction. The histograms 
were obtained from simulation data with diffusion coefficient 
D = 5x10"^ at volumes V = 10"* (red, unimodal), V = 2x10^ 
(purple, flat) and V = 10'^ (blue, bimodal). In each case the 
corresponding theoretical prediction of Eq. (|f 0|l is shown with 
a solid line. 



erful methods. This point was elucidated in Ohkubo et 
al flS^ , with the investigation of a simple one-dimensional 
model inspired by the TK reaction. 

We begin om- analysis of Eq. ([3]) for the TK reaction 
by making a change of variables which can be understood 
mathematically (as a real Fourier transform) or physi- 
cally (as corresponding to the total concentration, the z 
variable introduced by TK and two variables related to 
the Xi — X3 and X2 — X4 dynamics) . This is 

W = Xi + X2 + X-i + X4, Z ~ {xx + X3) — (X2 -f 2:4), 

U = X\ — X3, V — X2~ X4. (5) 

Applying the transformation, for the total concentration 
we find the closed equation w = D(4 — w). For the re- 
mainder of the analysis, we fix w to its fixed-point value 
of 4. For the variables z, u, and v we then find 
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where 4> = ^(z -I- 4)^/4 - , ^ = a/(z - 4)2/4 - v^, 
and the ( variables are independent Gaussian white noise. 

The dynamics on the 0(1) timescale are solvable. In 
fact, and tp defined above are conserved quantities of 
the system ([6]) with D set to zero. Solution trajectories 
are therefore confined to the closed curve given by the 
intersection of the surfaces defined by the values of (j) and 
ip, which are determined by initial conditions. Details of 
the full solution will be provided in a forthcoming paper 
[l4l] ; for the present discussion it is sufficient to point out 



that the trajectories are periodic, with the period for z 
being 



T 
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where K(- • • ) denotes the elliptic integral of the first 
kind. The period for u and v is double that of z. It 
is important to note that K{x) grows without bound 
as X — ^ 1, and thus Eq. ^ implies that the period of 
oscillation T diverges as either (/>—>■ or •(/;—>■ 0. In 
these limits, the trajectories of the deterministic dynam- 
ics deform into a homoclinic network linking the fixed 
point {u,v,z) = (0, 4, — 4) to (m, w, z) ~ (0, — 4, — 4) or 
(—4,0,4) to (4,0,4), respectively. This fact explains the 
presence of both fast oscillatory dynamics and metasta- 
bility in the same parameter range. 

We turn now to the study of the behavior of z on an 
0{1/D) timescale. From left to right, the terms in the 
equation for z in system (|6]) are responsible for the fast 
oscillation caused by interaction with u and u, the linear 
drag towards zero, and the noise. Since the oscillations 
occur on a timescale faster than the other two terms, we 
expect that a time average on a timescale r, such that 
T <^ T ^ 1/D, will not affect the drag and the noise 
substantially. To do the averaging, we coarse-grain time 
by intervals with length r in Eq. ^ and replace every 
term with its time average over that interval [15[. We 
write (• • ■ ) = Jt^'^ dt {■ ■ ■) for the time average and 
make use of the following assumptions: 
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0, (16 -z2); 
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These are justified on physical grounds: the first follows 
from the fact that the conserved quantities of the fast dy- 
namics are approximately constant on intervals of length 
r <C since the average of uv is a multiple of the 

average of z, and z has periodic trajectories if (j) and ip 
are fixed; in the second approximation, we are assuming 
that the strength of noise is not strongly affected by fast 
oscillations in z. 

The resulting so-called averaged equation for z is 



z = ~Dz 



Z-2)C(t) 



(9) 



This equation describes an interplay between the drag 
and the noise, and provides a complete picture of the 
dynamical transition first observed by TK. Physically, 
we may think of the system as gently relaxing to the 
origin, while being agitated by a noise term which van- 
ishes at the metastable states z = ±4. Depending on the 
strength of the noise (controlled by the parameter A), the 
system will either be attracted to zero by the linear drag, 
or forced to the boundaries by the noise. By varying A 
we transition between these dynamical regimes, an effect 
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which is most clearly demonstrated by calculation of the 
stationary distribution P{z; A). From Eq. (jO]), we find 



P(z; A) (16- r 



> A-l 



r(^ + A) 
^42A-ir(A)' 



(10) 



Our prediction is tested against the numerics in Fig. [2] 
This equation confirms the critical volume Vc = ^/D 
as the point of transition between a unimodal and bi- 
modal stationary distribution. We should point out that 
Eq. (|10|) is correct only up to first order in D; certain 
features of the simulation data (such as \z\ occasionally 
exceeding 4 due to variations in total concentration w) 
are not captured at this level of approximation. 

It is worth pausing a moment to discuss the relation of 
the noise- induced metastable states to the fast oscillatory 
dynamics discussed in the earlier analysis. For example, 
from the definition of V', we see that z = 4 can only be 
obtained when v = and "0 = 0, and thus we are in the 
regime in which the period of the oscillation is divergent. 
In this case one can expect fast periodicity to break down 
and the system to remain in a given metastable state for a 
random length of time, before being freed and proceeding 
along a trajectory close to the homoclinic orbit linking 
it to another. The lower right-hand plot of Fig. 1 shows 
this behavior. 

Beyond determining the stationary distribution of z, 
our methods may also be used to calculate various other 
uantities associated with the model. For example, in 
the fraction of time spent in the pair-dominated state 
(that is, Xi + X3 = or ^2 + X4 = 0), called the 'rate 
of residence', was measured from simulations and plotted 
as a function of A. The authors noted a puzzling shift in 
this quantity when adjusting for different cell volumes, 
which we are now able to explain. 

From Eq. ^ we can determine a straightforward pre- 
diction for the rate of residence by computing the fraction 
of time that z spends within 1/V of ±4. We integrate 
the stationary distribution to find: 



f4-l/V 

1-1 P{z; DV) dz 

-i+l/V 

-DV 
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V -\- higher order terms. 



Therefore, to properly compare different cell volumes and 
diffusion coefficients, one should hold DV ln{V) constant, 
rather than A. The fit between Eq. pT|) and data from 
simulations is shown in Fig. [31 

In this Rapid Communication we have examined the 
influence of noise on the link between structure and func- 
tion in a class of biochemical networks. The consistent 
formulation of the problem which we provide starts from 
the master equation and proceeds through a well-defined 
approximation scheme to an SDE which correctly cap- 
tures the behavior of the system. Although this equation 




DVln{V) 

FIG. 3. (Color online) Rate of residence of the pair-dominated 
state as a function of DVln{V). The circles show the result 
measured from simulations carried out with fixed D — 10~^ 
and varying V; for each data point a single simulation of 
duration tmax ~ 10^ was conducted and the fraction of time 
spent in the pair-dominated state measured. The solid line 
corresponds to the first-order prediction in Eq. (fTTj) . 



is not exactly solvable, we are able to proceed by identify- 
ing and exploiting a separation of timescales involved in 
the problem. This analytical process was demonstrated 
explicitly for the paradigmatic TK reaction, providing an 
understanding of the phenomenology of the model and 
yielding expressions for quantities of interest which are 
compared to the ones obtained numerically by TK. 

Since it is the discreteness of molecules which gives rise 
to the intrinsic noise experienced by reaction systems of 
this type, one might expect that such effects are only rel- 
evant in small systems, and can be neglected in general 
(indeed, this is a central assumption of any theory based 
on the study of macroscopic rate equations). In practice 
the situation is far more subtle; what matters more than 
the strength of the noise is how it interacts with other 
aspects of the model, such as the slow relaxation due to a 
small diffusion coefficient. As we have shown, this inter- 
action gives rise to metastability in the class of autocat- 
alytic reaction networks we investigate; moreover, it can 
be exploited mathematically to explain the dynamical 
transition observed in the TK reaction. A closely related 
noise effect has recently been observed in an ecological 
model [l6| . where it induces the spontaneous formation 
of species, and we expect that more surprising results of 
this type will come in the near future. 
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